This paper presents a method to identify gravitational arcs or more generally elongated struc- 
tures in a given image. The method is based on the computation of a local estimator of the 
elongation. The estimation of the local elongation proceed in two steps: first the local orienta- 
tion of the structure is computed, then in the next step, a rotation is performed, and the marginal 
distributions are used to compute the elongation. This procedure allows the computation of the 
local elongation at each point in the image. Then, using a threshold on the elongation map the 
elongated structures are identified and re-constructed using connectivity criteria. Finally a cata- 
log of elongated structures is produced, and the properties of each object are computed, allowing 
the selection of potential arc candidates. The final selection of the arc candidates is performed by 
visual inspection of multi-color images of a small number of objects. This method is a general 
tool that may be applied not only to gravitational arcs, but to all problems related to the mapping 
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and measurement of elongated structures, in an image, or a volume. 
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Abstract. 

1. Introduction 

The new large scale surveys offers deep and well sampled views of large regions of the sky. 
These views contains such a large number of objects that automated method to find a specific 
type of objects is of great interest. In particular, the new release of the CFHTLS survey (see 
\http://www. cflit. hawaii. edu/Science/CFHLS\ for a description of the project) presents, a large set 
of 1 square degree images, with a corresponding size of about 20000 x 20000 pixels. Among a 
large number of interesting objects, this survey contains a number of gravitational arcs, typically 
one per image. Considering the huge number of pixels, finding theses arcs or arclets is a good task 
for an automated software. The identification of potential arcs in the image is mostly a search for 
elongated objects. Thus an automated method designed to find arcs should first produce a catalog 
of elongated objects and select possible arcs candidates among them. In the end the procedure 
should be concluded by a visual inspection of a small number of candidates. This method is a 
new approach, former methods relied on filtering, for instance, Lerzer, Scherzer & Schindler, 
(2004 & 2005), used anisotropic filtering, while Starck, Donoho & Candes (2003) used wavelet 
Filtering. Arc detections using catalogs from the Sextractor software (Berlin & Arnouts 1996), 
were also performed by Horesh et al. 2005. 

2. The method 

In general arcs or arclets appear in the image as very elongated structures, almost flat along the 
elongated direction and nearly as narrow as the PSF along the other. Other structures in the image 
may have similar elongations, but along the smaller dimension they are rarely as narrow as the 
PSF. Some spiral arms or a random coincidence of small and faint structure, may be confused 
Send offprint requests to: C. Alard 
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with arcs, but most structures are not as narrow an arclet can be. Thus by selecting elongated and 
narrow structures in the image, we should build a sample which contains a significant fraction of 
arcs. Once this sample is constructed the next step will be to characterize the arcs in the sample, 
with respects to the other contaminants, which can be real structures in the image, but also defects 
due to bright stars. The paper is organized as follows: in the first part, an estimator of the local 
elongation at small scale is introduced. The second part analyze the effect of the noise in the 
image on the estimation of the local elongation. And finally, some illustration of the practical 
implementation of the method will be given. 



2.1. Reconstruction of the narrow elongated structures in the image 
2.1 .1 . Estimating the local elongation at small scale 

At each point in the image defined by its coordinates (xq, Jo) we wish to estimate the local elon- 
gation at the scale of the PSF. The elongation is computed by analyzing the pixels in a window 
of size about a few times the eflFective size of the PSF around the point of interest. The first point 
in this analysis is to estimate the local orientation of the structure, by computing its second order 
moments. Then using the orientation defined by the moments, the local axis ratio is estimated us- 
ing the 2 marginal distributions /y, and, Ix in the proper axis. To be more specific, the orientation 
is defined by the second orders moments: 

a = J lixo + xuyo+ yi) xl dxidyi 
b = ^ Kxo + xuyo+yi) y\ dxidyi 

c = J lixo + xi,yo xiyi dxidyi 

The coordinates (xi,yi), are local coordinates with origin at the image coordinates {xo,yo). The 
integration are performed in the interval -M < xi < M, and -M <y\ <M. The direction of the 
main axis is then defined by the angle: 

tan(0) = \ , , = (1) 

-1/2^-1- l/2a- 1/2 ib^-2ab + a^ + Ac^ 

We perform the rotation with angle 6 and center (xq.Jo). in order to move from the original 

local coordinates {x\,yi), to the proper coordinates {x,y). In proper coordinates the direction of 

maximum elongation is along the x axis. We then define The marginal distribution: 



Ix{x) = J 

iy(y) = J 



I(xo + x,yo+y)dy 



I(xo + x,yo + y)dx 

At each position in the image (xo,>'o), we define the following measurement of the local elonga- 
tion: 

1 /f(0) 

Q{xo, yo) = TT7 pTfV (2) 

2M sup LIx(x)]i-M<x<M] 
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each direction, then locally, /(xq + x,yo+y) = fix) g(y) and we infer that: 

1 r f(x) dx 

2M sup [/(x)][_^<^<^ 

In this calculation it was assumed that the profile along the small axis g(y) is normalized, thus 
that Jg(y)dy=l. 

2.1 .2. Properties of the estimator 

The first case of interest is an arclet with flat profile along the tangential direction. Then, Ix is 
constant, and we have simply: 

Q(x,y) = g(0) 

Assuming that the profile g(y) may be written as a generic function of the width of the profile a, 
g(y) - \ G with G{u) symmetrical in u and thus maximum at w = 0, we infer immediately 
that the estimator Q{x, y) will be maximum at the center of the arc with a value Q{x, y) = \ G(0), 
which means clearly that Q{x, y) wiU reach its maximum value at the center of the narrowest 
arcs. We may turn now to the general case, when the distribution Ix isn't flat. It is easy to find an 
upper bound for the integral on f{x) in Q{x,y) by using the upper bound on Ix, since sup/x = 
sup [/(x)], we infer that: 

Q{x,y) < g{Q) 

Thus this estimator wiU be always maximum on flat and narrow elongated structures. Any varia- 
tion along the tangential direction may be considered as the consequence of a flnite scale length 
along this direction. In particular if the scale length of the variation is smaller than M, the esti- 
mator Q{x, y) wiU give directly an estimation of the axis ratio. Assuming: 



And b < M, then 



Q(x,y)=kbg(0)^k-G(0) 
a 



When b is larger than M the estimator Q{x, y) is no longer the axis ratio. Actually the value of the 
main axis b saturate at the value of the box size M. The Umiting value M for b is reached when 
the distribution becomes flat or equivalently when the efl'ective scale-length becomes infinite. 



2.1 .3. Effect of the noise on the estimator 

We may now evaluate the fluctuation of the estimator Qix,y) due to the noise in the image. 
The noise on the estimation of the local axis ratio may be decomposed in 2 parts. First, given 
an orientation of the local axis, there is some amount of noise due to the intrinsic noise on the 
marginal distribution. This part wiU be computed first. But there is also another source of noise 
(a priori un-correlated), which is due to the noise on the orientation 6 of the local structure. 
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Assuming a variance o-y for ly, and (Tx for sup [Ix], the variance on Q may be expressed as: 



oy (Tx ^ 
\'y ^x ^ 



The variance cry of the marginal distribution lyiy) scales like yw+i ^^'-^ respect to the variance 

on the total flux in the box. The signal will scale like the integral on the profile marginal 

distribution giy) in a bin of with 1 pixel normalized by the total integral on g{y). Since the arcs 

have about the width of PSF, and that the PSF may be modeled by a gaussian of width 4 pixels, 

we find that Iy{x) at the maximum of the arc will be about 1 /4 of the total signal S . Then: 

(Ty ^ 16 

If ~ 2M+ 1 52 

We turn now to the estimation of <tx. Considering that the distribution probability of the noise 
has a partition function 0{u), and a probability distribution <f>{u), the probability distribution of 
the maximum of the distribution will be: 

Pm(m) = (2 M + 1) 4>iu) ^{uf^ 

Assuming a gaussian probability distribution for the noise, we find that: 

Pm(u) = K exp(-u^/a^) (1 + erf(u/o-)) 

This distribution may be approximately modeled by a gaussian distribution with an ofl'set uqs: 

Px{u) = A exp(-(M - Mo)^/cr). for M ^ 1.5 seeing ^ 7 pixels, we find that a ^ 1/4 cr. And 

finally we infer the signal to noise ratio the variance o"x: 

Ox ^ i2M+l)N^ 
ll ^ 452 

And eventually, taking M = 7, for a well sample PSF, the signal to noise ratio for the estimator 
Q reads: 

Q Q 



Nq ylo^ 

Which means that even for low signal to noise, ^ ^ 10, the noise will be only 20 % of the signal. 
This is small compared to the usual dynamics on Q itself, which for our choice of parameter 
varies from about 1 to 3. 



2.1 .4. Noise on the estimation of the local orientation 

An important problem is to estimate the error due to the noise image that we make in the estima- 
tion of 6. Most of the noise on the arclets is due to the fluctuation of the background, thus it will 
be assumed the amplitude of the noise is nearly constant. Lets call a the variance of the pixels 
due to the background noise, the relevant noise on the moments is then: 

<a^>^Yj< ^Ij ^(T r / dxdy = a 2/5 

Similarly: 

<b^>=o- 2/5M^ 
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And: 

< >= 0-4/9 

In order to obtain the variance on it is necessary to differentiate equation O with respect to 
(a,b,c), and to add the respective variances. The final result may be expressed using the two 
proper second order momentum (a, and /3) of the system. Since: 

a = sin(6')^ Q3 - a) + a 

b^sinief {a-l3)+l3 

sin(20) , „ 
c = — - — {a-p) 



, cos {ley 2 M'' 



A useful quantity is the average variance. 

M^cr 

<02>- 1/5. 



Considering that the total noise in a window of size M is: 

= cr(2M+ \ f 

And that the proper second order momentum may be rewritten, as a function of the proper size 
Lo, L\ , and of the total flux in the window, S : 

It is possible to estimate the variance on 9: 

< 0^ > T- 

For a flat distribution like an arclet the size Lq along the main axis is close to the size of the 
window M, while the other dimension is small and may be neglected. In this case the former 
expression reduces to: 

<02>.1/2O — 

Which means that the error on 6 depends only on the signal to noise. For small signal to noise, 
typically ^ - 10 the error on 6 is only of about: 



N(e) = V< 02 > ^ Vl/200 radians ^ 1" 



The error on 6 has an effect on the estimation of the local elongation. The arclet will be repre- 
sented by a quadratic polynomial in proper coordinates: 

I(x,y) - F(ao + ai y^) 
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Assuming that the arclet is mis-aligned by an angle dO with its proper axis, the intensity of the 
arclet now reads: 

I = F ([ao + {fli - ao)d(p-^ + [ao + (ao - ai)<i0^] + 2 {oq - ai)de ry) 

In particular we are interested with the effect of this rotation on the scale of the distribution along 
its small axis. The distribution along a particular axis may be estimated by using the reduced 
one dimensional marginal distribution along this axis. Thus the effect of the rotation may be 
estimated by calculating the effect on the marginal distribution along the y axis: 

ly = J" I{x,y)dx 

Which may be re- written by changing the integration variable: 

h = J F{au^ +liy^)du 

WithyS: 

yS = ai(l-i-(l-ai/ao) dO^) 

By introducing the scale length Lq = 1 / Voo, and L\ - Ij -\/ai, and the new scale length of the 
marginal distribution, Li - 1/ we obtain: 



U=LA\ + - 



-1 + =^ 



d6^ 



As estimated in the former section, dO^ ^ and in in our choice of parameters, ^ ^ 3 , we 
find that the relative difference between L\ and L\ is only about 2 %. The effect of the rotation 
on the main axis smaller again, and thus the error on Q has a negligible effect on the estimation 
of the axis ratio. 



2.2. Effect of local curvature 

In this paper it was assumed that arclets or elongated structures may be de-composed in small 
rectilinear bits at the scale 2M, which is basically a few times the size of the PSF. This is an 
excellent approximation for larger arcs, and the estimator Qix,y) is designed for this specific 
task. But it is interesting to study the behavior of the estimator for smaller arcs. As usual we 
assume that without curvature at local scale the profile can be de-composed along each direction: 
I{xo +x,yo+y) = fix) g(y). Assuming that structure is curved locally and may be approximated 
by a circle of radius R and that the tangential direction is along the x axis, the curvature introduce 
a small shift 6y of the profile along the y direction. To the lowest order this shift is: 

The shift 6y has a negligible effect on Ix, thus the effect on Q{x, y) will be related to the variation 
due to SY on ly only. With curvature, ly may be re-written: 

hiy) = j fix) giy + 6y)dx - g(y) j f{x)dx + ^ J f{x) ^dx + ^ J fix) ^dx 
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Assuming that the profile along the y direction, g{y) is symmetrical and that the derivative is con- 
tinuous, we have ^ = 0. For an arclet, the profile along the x direction /(x) may be approximated 
by a constant /o in the relevant interval. Then /^(O) reads: 



iY = 2Mgm\ + 



40/?2 



\_<fg 
gcfy 



y=0. 



We have now to evaluate the second derivative of the functional g near y = 0. Expanding g to 
order 2 in j, we have: 



g(y)-A\i + - 



and subsequently: 



/y^2Mg(0) 1 + 



1 /M\2/M\2\ 



Numerically, the perturbation introduced by curvature may be evaluated as follows: M is about 
1 .5 times the seeing, S , and since a = 2 VC log(2))5' , we have, M ^ 2.55 . And finally we have: 



A.. 1.95(1) 



Basically, this means that the error is still acceptable (20 %, about the noise level for faint ob- 
jects), when, R, is about 3 times the seeing. 



3. Practical implementation. 

3.1. Filtering 

Finding arcs is equivalent to look for structures that have nearly the size of the PSF along one 
dimension. There is no scale length smaller than the PSF in the image, except artifacts due to the 
noise. Thus any filter that damps the lower frequencies while damping somewhat the noise will 
improve the detectability of the arcs. Such filters are basically band-pass filter, centered at the 
PSF scale. An example of such filter is a Mexican-hat filter: 

M{x,y) = exp -0.5 exp ^ (3) 

The value of the constant b should estimated so that the filter width match the size of the PSF. 
It is clear that this 2-dimensional filter, will improve the detectability of structures at about the 
size of the PSF. But the real point of interest is the efl'ect of the filter on a one-dimensional 
structure like an arclet. To estimate the efl'ect of the filter on the arc, we may take advantage of the 
symmetry of the filter by making a local rotation around the filter center such that the X-axis will 
coincide with the direction tangential to the arc. Assuming that the arc has a a flat profile along 
the tangential direction, and a gaussian profile along the narrow direction: I{x) = exp{-x^/a^), 
we may calculate the scalar-product with the mexican-hat defined in Eq ( 1 ) . The amplitude of this 
scalar-product will be proportional to the amplitude of the arc after filtering. The scalar-product 
reads: 

P(/7)=;r77cr2 — 
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Number of pixels 


mean value 


pix nb/radius 


radius 


remarks 


103 


2.19 


3.8 


26.8 


Central arclet 


59 


1.93 


1.0 


56.4 




53 


1.81 


0.9 


58.5 




46 


2.15 


0.9 


50.2 




45 


1.59 


0.8 


50.7 





Table 1. paramaters of the structures found in the Q(x,y) image (Fig.|2j- 



Number of pixels 


mean value 


pix nb/radius 


radius 


remarks 


81 


2.13 


1.90 


42.7 




79 


1.97 


1.28 


61.6 




65 


1.88 


1.16 


55.8 




58 


1.61 


0.91 


63.5 




45 


1.91 


3.14 


14.3 


Small arclet 



Table 2. paramaters of the structures found in the Q(x, y) image. 



The functional P(t]) peaks at 77 ^ 0.838, and decrease quickly for larger values of 77. Thus it is 
clear that for maximum efficiency, given a scale a for the PSF, the scale cr of the filter must be: 
o- ^ 1.19fl. 

3.2. Using the estimator Qix,y) to find elongated structures. 

As explained in the former sections, at each point in the image the direction 6 of the local structure 
is computed using the second order moments, a rotation is performed in order to be aligned 
with the direction of the local structure and finally the estimator Q(x, y) is computed using Eq. 
(|2j. To illustrate the method, an image representing the value of the estimator Q{x,y) at each 
point of the image presented in Fig. Q is presented in Fig. Local elongation of about 3 are 
reached within the arclet at the center of the image, and no-where else in the image. Actually 
if larger images were taken, some objects other than arclets, may also presents large value of 
the local elongation. For instance, edge-on spirals. Then the problem is to separate the arclets 
from the rest of the sample. This can be done mostly by computing the properties of the local 
structures as detected in Fig. (|2}. The structures in the Q{x,y) image may be reconstructed using 
the ink-blot technique. A structure is defined as a set of pixels larger than some value, and the 
pixel must be connected each other. This way all the structures in the Q(x,y) image may be 
reconstructed,and their properties: number of pixels, length, curvature,. ..ect, may be measured. It 
is generally possible to separate most of the arclets form other structures using these parameters. 
For this particular task it may be also useful to measure some of the parameters of the structures 
in the original image, like for instance their color or brightness. Table (1) gives the 5 largest 
structures found in the image presented in Fig. (|2j- Note that the actual field that was used is 
much larger than in Fig. (|2}, the effective field was of 1000x1000 pixels for Table 1. 
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100 200 300 

X (pixels) 

Fig. 1. Contour plot of a CFHT image. In this figure an arclet is visible at the center of the frame, 
near a bright source. This image will be used to illustrate the computation of the local estimator 
Qix,y)- (sec Fig. 2). 

3.3. Computing time 

Considering the size of the image, the computation of the estimator Qix,y) as to be fast. The 
filtering and the computation of the second order moments can be greatly optimized using well 
known techniques, like the decomposability of the filters, or partial buffering of the calculations. 
However, it is more difficult to speed-up the calculation of the marginal distributions and of the 
estimator itself, since the orientation is not known in advance. But it is clear that the computation 
is not meaningful in all points of the image, first it is possible to introduce a threshold in the 
filtered image to perform the computation. A further gain of computing time may be operated 
by calculating the estimator only near uni-dimensional local maxima. Uni-dimensional local- 
maximas are maximas along a particular direction, like for instance the kind of maxima that will 
be found on a direction perpendicular to an arclet. Since a priori the orientation of the arclet is 
not known, the maximas wiU be searched both along the x and y axis. Each local maxima along 
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Fig. 2. The contour plots in this image represents the local values of the estimator Q{x, y) at 
each point of the image presented in Fig. 1. The contours were plotted using a linear scale, and 
the highest contours (which are located in the arclet) are about 3 times higher than the lowest 
contours. 

the X or y axis will be defined as a point where the estimator Q,{x,y) has to be calculated. This 
procedure speeds-up the computation of Q.{x,y) by about a factor of 10 in a typical CFHTLS 
image. 
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Fig. 3. A small and faint arclet is located near the center of the image. Note that this small arc is 
very close to the deflector, and that there is some blending of the light with the lens. 
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Fig. 4. When estimating the local elongation Q{x, y) the arclet appears now clearly, and stands 
out in the image. 



